COMPUTATIONAL ASPECTS OF THE PREDICTION OF MULTIDIMENSIONAL 
TRANSONIC FLOWS IN TURBOMACHINERY 

By David A. Oliver and Panagiotis Sparis 
Massachusetts Institute of Technology 

INTRODUCTION 

The analytical prediction and description of transonic flow in turbomachinery is 
complicated by three fundamental effects: (1) The fluid equations describing the transonic 
regime are inherently nonlinear, (2) shock waves may be present in the flow, and (3) turbo- 
machine blading is geometrically complex, possessing large amounts of curvature, stagger, 
and twist. Simple analytically separable solutions are therefore not readily obtainable. 
(The complex geometry of a t 3 rpical transonic compressor rotor is shown in fig. 1.) 

Because of these analytical difficulties, a computational approach to the prediction and 
design of transonic turbomachine flows is strongly warranted. 

In the present work, a three-dimensional computation procedure for the study of 
transonic turbomachine fluid mechanics is described. The fluid differential equations 
and corresponding difference operators are presented, the boundary conditions for com- 
plex blade shapes are described, and the computational implementation and mapping pro- 
cedures are developed. Illustrative results of a typical unthrottled transonic rotor are • 
also presented. 

’ FLUID EQUATIONS AND DIFFERENCE OPERATORS 

The densities of mass p, momenta mj, and energy e defined by the fluid state 
vector U(xi,t) are governed by the fluid conservation laws in cylindrical coordinates 
^xj = r, X 2 = 0 , and X 3 = z) and time t: 



In a coordinate system fixed to a turbomachine blade rotating at angular velocity 12, 
the state vector U, and the flux vectors Fj^(U) and K(U) are 

U= [p,pur,pu0,puz,e] (2) 
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In equations (3), the pressure p appears which is expressed in terms of the state vec- 
tor U through the equation of state. 

A time-e:iq)licit difference operator S approximating the fluid equations (eq. (1)) 
with second-order accura.cy coupled with a local stabilizing operator D is used to , 
advance the fluid state from time level n to level n -f 1 over the interval 6t as 
follows: 


= (S + D)U" 


(4) 


Here, following MacCormack (ref. 1), S is formed in two steps as 

U* = - ^ AjJ:i(U) + 6«(uj] 

s« * 2 - Aj-F,®*) * 


(5) 


The operators and Aj are the forward and backward dlHerence operators in 
the coordinate directions Xj: 

A"''f(x) = f(x + 6x) - f(x) 


A"f(x) = f(x) .- f (x - 6x) 


( 6 ) 
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K the difference operator (eq. (5)) is applied successively in split steps (ref. 2) for 
each coordinate direction, the numeric^ stability conditions are 


6t ^ Min/— (7) 
i |[ui|+.c[ 

At stagnation points or sonic points, the linearized version of the difference opera- 
tor (eq. (5)) is neutrally stable independently of the choice of 6t/6xj. For a genuinely 
linear difference operator this occurrence is of no consequence and stable solutions of 
the difference equations can be achieved. In the nonlinear case, however, the tine sta- 
bility of the operator at the neutral point will be determined by the higher order nonlinear 
terms, and these terms will destabilize the difference operator. This numerical instabil- 
ity is an inherently numerical instability induced by higher order terms in the tmncation 
error. Hence the nonlinear difference operator may be stabilized by the introduction of 
an artificial dissipation term of the order of the truncation error of the difference opera- 
tor. In the case of the difference operator (eq. (5)), this term must be of third order in 
6t,fiXi. 


The global damping or stabilizing operator Dq may be summarized in the arche- 
typal form 

> 


DoU 


6t a~[qa[ + a^qaI 

6xi 2 


(8) 


where Q(U) is a matrix diffusion coefficient of order 6t,6xi and k is an order unity 
nondimensional constant. A Taylor's series expansion of this operator shows that it is 
the difference expression for the continuous diffusion operator 


° 1 6Xi3Xi^axi 


(9) 


The neutral stability condition of the linearized operator occurs at sonic points and stagna- 
tion points. While a sonic point exists in the interior of a shock wave, the nonlinear insta- 
bility could be just as important in smooth isentropic regions of transonic flow passing 
continuously through the sonic point as in regions where shocks are present. Calculations 
performed with the operator (eq. (5)) have confirmed this conjecture. For subsonic flows 
up into the high subsonic regime, the undamped operators have been demonstrated to be 
stable. However, where the Mach number was increased into the transonic regime, numer- 
ical instabilities occurred which could be eliminated with the use of the damping operator. 

On the basis of the observations, the nonlinear instability should be confined to those 
regions of the flow where the eigenvalues of the amplification matrix are unity, i.e., at 
sonic and stagnation points. The damping operator D should then be structured such 
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that it'becdmes significant only near the sonic and stagnation points and not operative in 
other portions of the flow. A damping operator with such characteristics is 

D = f(M)Do ° (id)" 

<r; > ' f.- J • 

where f(M) is a distribution function which depends upon the local Mach number such 
that f(M) = 1 for M = 0 ^or M = 1, but f(M) « 1 for M 0 and.l- A useful function 
with such properties is the Lorentz Line shape function 


f(M) = 

1 + 


M - 1 
AM 



( 11 ) 


where the paranieter AM represents the effective, width in Mach number of. the distriT 
button function which peaks at M = 0 and 1. ■ The. use of the local damping operator D. 
in place of the global operator Dq can offer significant improvement in the resolution, 
of the flow while maintaining stability of the difference operator. 

. The utility of the local damping operator is shown in figure 2 where the supercriti 
cal transonic flow of a y = lA gas over a right circular cylinder has ..been computed. 
This calculation was performed with a minimal number of mesh points (20) distributed 
over the surface of the half- cylinder to test the utility of the local damping operator. 
Si^ificant improvement of the flow -field resolution, including the shock wave, results. 


BOUNDARY CONDITIONS ON BLADE SURFACES 


The appropriate boundary condition on an. impenetrable blade surface for an invis- 
cid flow is the single condition 


- 3 un = 0 (12) 

where Uj^ is the velocity normal to the blade surface in blade coordinates. This bound- 
ary, condition is not readily implemented in the finite-difference procedures described in 
the previous section because the full fluid state U is required at each point includii^ 
the boundary points. If, as in figure 3, the boundary points are treated as interior points 
to which the difference equation (eq. (5)) is to be applied, then the application of the bound- 
ary condition, which consists of the determination of the state vector U at the auxiliary 
point, must be such that only the single condition (eq. (12)) is imposed'at the blade surface. 
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If 2 is the surface of a sufficiently smooth three-dimensional body, then at each 
point M on the surface 2 a triply orthogonal curvilinear coordinate system may be 
defined consisting of the local normal ^ to the surface 2 at the point M and two 
cuzres of the surface 2(77 and |) that are normal to each other at the point M. If 
dn, dr, and ds are the differential arc lei^ths along the axes rj, and respec- 
tively, then 

ds=hsdri 


dr = hrdq\ 
dn = hn d^ 


(13) 


The equations of motion may be expressed in the curvilinear coordinate system. 
q, and C* By introducing the boundary condition U{) = 0 and the radii of curvature Rs 
and Rj of the surface 2 in the direction of the axes | and 77 , the equations of 
motion may be written in the following form, correct only for points on the surface '.';2 
(ref. 3): 1 


8t 
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8t hs 8| hr 877 hshr 877 hshj 8^ phs 8^ 


(14) 
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(17) 


Equations (14) to (17) have been obtained with the isentropic assumption in which a^ is 
the square of the speed of sound. By eliminating the space and time derivatives of the 
density between equations (14) and (17) and replacing the general curvilinear coordinates 
|,_77, and ^ with the local s, t, and n, the following condition on the normal deriva- 
tives of the normal velocity is obtained (ref. 3): 
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By center differencing the normal first and second derivatives of the normal veloc- 
ity, an approximation is obtained to u“ on the auxiliary point correct to order (6n)^ if 
the right-hand side of equation (18) can be approximated with an error no greater than 
0(6n), Thus, on the surface of the body, equations (14), (15), and (16) have to be solved 
while a boundary condition is applied for u^ using equation (18). Note also that the 
radii of curvature of the surface Rg and R.^ are required for this accurate determina- 
tion of uj^. The second normal derivative of the normal velocity is rigorously required 

for a second-order-accurate boundary, condition, the expression for uj^ beii^ 

. * . , • - ■ .... - 


However, for mildly curvii^ shapes with. 1/R of order 6n,, a simple reflection of u„ 
will yield second-order accuracy without the complication of introducing the second nor- 
mal derivative given by equation (18). This simplification is used in the results to be 
illustrated in the subsequent sections. 

UPSTREAM, DOWNSTREAM BOUNDARY CONDITIONS 

For a steady-state one -dimensional duct flow, one is not free to specify the down- 
stream boundary conditions if the upstream conditions are fully specified. This is not so 
in a transonic multidimensional duct flow. In this case the upstream conditions may be 
set. However, one can still vary a single downstream variable such, as the pressure and 
achieve different steady-state solutions. This degree of freedom on downstream pressure 
arises because the oblique shock waves present on the blades in the transonic regime are 
free to move and alter their strength in response to the different downstream pressure 
conditions. This range of freedom on downstream pressure is limited. It ceases, for;;. 
example, if the downstream pressure is set high enough so that the shocks are blown for-; 
ward out of the cascade. ‘ . - • , , . 


In the actual calculation, boundary conditions were set so that the mean flow at the 
inlet plane to the duct was held fixed. Waves generated by the rotor which ihbved' 
upstream into the inlet plane were then allowed to escape. This escape condition was 
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formulated as an axially one-dimensional characteristics construction at the inlet plane 
Of the computational domain. At the downstream exit plane the pressure was held fixed 
and the remaining flow variables forced to take a zero axial gradient condition. This 
condition allows the mean flow velocity at the exit to adjust itself to the correct mass 
flow; however, it distorts the structure imposed by the rotor loc^ly in the vicinity of the 
^t plane. This distortion is a consequence of the condition of uniformity of flow in the 
axial direction rather than in the streamline direction (which is helical rather than axial). 


COMPUTATIONAL IMPLEMENTATION 

The geometry of the flow field for an illustrative transonic rotor calculation is 
shown in figures 4 and 5. On the conical spinner are attached N = 23 blades with an 
average hub to tip ratio of 0.6. Thus the flow field of a blade element is bounded by the 
machine outer and inner casings in the radial direction and by two surfaces 0j(r,z) and 
02 (r,z) in the angular direction such that 


. V 02(r,z) - 0i(r,z) = A0 , ... 

A0 = (20) 

L.. 

Let rjj(z) and r,j,(z) be the equations of the inner casing (hub) and the outer 
casing (tip) surfaces. If 0 = 6q{t,z) and 0 = 0p(r,z) are the equations of the blade 
suction and pressure surfaces, respectively, then the computational domain boundary sur- 
faces 0j and 02 may be constructed by making 0i = 0s and ^ 2 = 9 ^ in the blade 
region and then extending these surfaces upstream and downstream as ruled surfaces 
parallel to the machine axis as shown in figure 4. 

The complex geometry formed by the extended blade and casing surfaces may be 
handled computationally by carrying out the computational work in a computational domain 
xj = (r',0’,z*) obtained as a mapping of the physical domain xj = (r,0,z). The complex 
tiirbomachine surfaces should map into planar surfaces in the computational domain. 
Mapping functions selected for this purpose are 


r - rg(z) 


^T 


(21) 


0 ' = 


0 - 0l(r,z) 
02(r,z) - 0i(r,z) 


A0 


( 22 ) 
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(23) 


2 = Z 


8X1 


The Jacobian derivatives of this transformation 




8x. 




appear in the conservation 


8F. 8F. 

I-.— ^ 


laws (eq. (1)) so that — - is replaced by g. . — =-. Similarly, the difference ojperator 

8Xj . 8x1 

(eq. (5)) in the computation domain becomes 


U* = - 


r(r’,z’) 


gij ^111(11) + StK® 
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(24) 


In addition, the blade boundary conditions which require the normal derivatives in 
the physical domain must be expressed in terms of computational domain coordinates. 

If v^ . denotes the three direction cosines of the blade surface normal, then the first nor> 
mal derivative expressed in computational domain coordinates Xj is 


8 8 

8n " 8x. 


(25) 


and the second normal derivative is 



*i*^k®ij ®kl 


8 


( 26 ) 


The derivative forms (eqs. (25) and (26)) and similar forms for the s and t deriva 
tives replace those of equations (14) to (18) allowing finite differences to be taken in 
terms of computational domain coordinates xl. 

INIlivL RESULTS FOR A TRANSONIC ROTOR 


Some initial results obtained with the foregoing method for the single-stage tran- 
sonic rotor shown in figure 1 are now presented. This rotor operates with a tip M&ch 
number Mij. of 1.2 and an average axial Mach number of 0.5. The calculation 
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illustrated here is for an open throttle situation in which the static pressure behind the 

k i. 

rotor remains at a low value relative to the full load design values of the rotor. 

. The calculation was performed on an IBM 370/168 system, i. Although this machine 
has a large primary memory capacity, economic considerations necessitated the use of 
the secondary memory units for the storage of the and arrays. Continuous 

input-output operations between the core and the secondary memory units were required 
for the calculation. . 

The computational domain shown in figures 4 and 5 was discretized with 67 points 
in the axial direction, 30 of which are in the^region of the blade. The radial direction 
was discretized with 12 points, the angular Erection with 10 points. This discretization 
was selected as being the minimum number of points which would provide a representa- 
tive although not detailed resolution of the flow field. Considerations of economy dictated 
the use of such a coarse 'mesh in a developmental calculation such as the one described 
here. Future calculations may be carried out with expanded, mesh densities. 

Before going into the details of the obtained results, it might be useful for the reader 
to familiarize himself with the nature of the coordinate transformation and especially with 
the shape of. the r' = Constant and 6' = Constant surfaces in the physical space. The 
r' = Constant surfaces are cylindrical surfaces as indicated in figure 5, ranging between 
the hub and the tip. Actually the tip and the hub surfaces belong to the -r' = Constant 
family of surfaces. The 0',= Constant /surfaces are more complicated and, as it can be 
seen in figure 4, they have no degree of symmetry. However, it is quite apparent that 
tiiese surfaces are more geometrically related to the blade shape than the 9 = Constant 
planes. Thus, some of the plots of the field properties have been made in the r', .0', 
and z' coordinate system rather than the r, 0, and z system to improve their 
clarity. . 

In figure 6, the Mach number distribution is plotted in r and z coordinates along 
,the span on the 0' = Constant surface that passes at the suction surface of the blade 
(0’ = A0). 

There is a strong shock wave originating from the trailing edge of the blade; on 
the other hand, there is hardly any trace of the leading-edge oblique shock. The flow 
upstream of the blade varies almost linearly from subsonic close to the hub to super- 
sonic at the tip, as expected, due to the solid body rotation of the flow in the rotary frame. 

Because of the wide open throttle condition, the flow accelerates near the hub to a 
higher Mach number (approximately 1.9 compared to 1.6 at the tip); This causes the 
shock wave at the trailing edge to be strongest near the hub. The flow at the supersonic 
tip is relatively smooth, the trailing-edge shock decelerating the flow* from a Mach num- 
ber of 1.6 to 1.3. 
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In figure 7, a similar plot of the Mach number distribution on the pressure side of 
the blade is illustrated. The effect of the hub geometry becomes more apparent now 
since the flow on the pressure side is accelerating. The maximum Mach number on the 
pressure side is 1.55 anil occurs at the tip 'j^ectiori. A shock is present also on the pres- 
sure side; however, its strer^th is diminishing towards the hub, where the flow becomes 
nearly sonic. The pressure side' shock appears upstrea;m of the trailing edge and, as will 
become apparent in the next figures, it extends to the trailing edge of the suction side of 
the adjacent blade. 

In figures 8 and 9, the intersections of the sonic surface with the pressure and the 
suction side have been traced in r and z coordinates. In figure 8, corresponding to 
the pressure side (0 = 0), the presence of the trailing-edge shock is rather apparent. 

This shock is normal to the hub surface at the pressure side, changing to an oblique 
towards the suction side (fig. 9), In figures 10, 11, 12, the Mach number contours 

have been plotted in 0 and z coordinates on the three r' = Constant cylinders cor- 
responding to the tip, mid, and hub sections of the blade. The location and the strength 
of the trailing- et^e shock can now be easily established. 

Figures 13 and 14 show the pressure coefficient 
r -P ■ Poo " ' 

^P ~2 

p a„“ 

^00 00 

(where ( )oo indicates the properties far upstream) on the suction and pressure sides at 
typical blade sections near the tip and the hub, respectively. The location of the trailing- 
ec^e shock can be easily established on both blade surfaces. 

In figure 13, the leading-edge shock can be clearly identified. This shock appears 
to be stronger on the suction side of the blade. At the pressure side near the trailing 
edge, there is evidence of an expansion fan that matches the local static pressure at the 
two sides of the blade. With respect to the blade's lift distribution along the span, the tip 
section produces little work at the open throttle condition; on the other hand, the hub sec- 
tions are working very hard, but in vain, in view of the strong shock at the trailing edge 
that dissipates much of the stagnation pressure rise produced at the hub. 

In figure 15 is plotted the distribution of the crossflow velocity Uj. on the 
0 = Constant surface that passes at the blade's suction surface (0 = A0). The behavior 
of Uj. is generally determined by the shape of the hub and tip; however, in the region 
around the trailing edge of the blade, die appearing shock has a significant effect. Due to 
the acceleration of the flow, the pressure upstream of the shock is lower at the hub than 
at the tip; thus, a crossflow is developed from the tip to the hub. On the contrary, down- 
stream of the shock, the crossflow reverses direction because now the hub pressure is 
higher since the shock is stronger at the hub. 

576 



Figures 13 and 14 also clearly show the acceleration of the flow, even in the super- 
sonic region near the blade tips, due to the influence of the subsonic zone near the hub. 
This is a crucial effect in such a genuinely mixed flow; for even though axially propagating 
signals cannot influence the supersonic zone, it is possible for disturbances to feed back 
upstream through the subsonic zone and then radially impact the upstream supersonic 
flow at the tip. 
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Figure 1.- Single stage transonic compressor rotor designed 
with average hub-tip radius ratio of 0.6; axial Mach num- 
ber, 0.5; tip Mach number, 1.2. 



Figure 2.- Test of local damping operator D(M) for 
transonic flow of y = 1.4 gas over right circular 
cylinder. Twenty points distributed over 6. 
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/VIESH POINT INDEX 


Figure 5,- Computational domain (section view) along axis of machine 
illustrating r' = Constant surfaces. 
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BLADE SURFACE 


Figure 6.- Mach number distribution over blade surface (suction side). 
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Figure 8.- Sonic surface intersections with 0 = surface (pressure side). 
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Figure 13.- Pressure coefficient through passage and over blades at tip. 

(1 in.; = 2.54 cm.) 



Figure 14.- Pressure coefficient through passage and over blades at hub. 

(1 in. = 2.54 cm.) 
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